The tunneling density-of-states of interacting massless Dirac fermions 
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We calculate the tunneling density-of-states (DOS) of a disorder-free two-dimensional interacting 
electron system with a massless-Dirac band Hamiltonian. The DOS exhibits two main features: i) 
linear growth at large energies with a slope that is suppressed by quasiparticle velocity enhancement, 
and ii) a rich structure of plasmaron peaks which appear at negative bias voltages in an n-doped 
sample and at positive bias voltages in a p-doped sample. We predict that the DOS at the Dirac 
point is non-zero even in the absence of disorder because of electron-electron interactions, and that 
it is then accurately proportional to the Fermi energy. The finite background DOS observed at the 
Dirac point of graphene sheets and topological insulator surfaces can therefore be an interaction 
effect rather than a disorder effect. 
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Introduction: A l^-shaped density-of-states (DOS) curve, 
vq{E) oc \E\/v'^ where E is energy measured from 
the Dirac point, is a distinctive feature of graphene 
and of topologial insulator (TI) surface states. Both 
two-dimensional electron systems are described by free- 
particle massless Dirac fermion (MDF) Hamiltonians: 
ho{k) = hvcT • fc [Hig. Here cr = (fi^, cr«) is a 2D vec- 
tor of Pauli matrices and v is the MDF velocity, which 
is roughly three-hundred times smaller than the speed of 
light in the graphene sheet case [3 [3] , and six- hundred 
times smaller [2^ in the case of the currently-studied 
chalcogenide TPs. 

The DOS of electronic systems has often been mea- 
sured using tunneling, and in recent decades often using 
scanning tunneling microscopy (STM) which probes lo- 
cally to mitigate the influence of inhomogeneities and 
some types of disorder. In STM [3 the DOS is propor- 
tional to the differential conductivity dlfip/ dVhias, where 
/tip and Vbias are respectively the tunneling current and 
the bias voltage between the STM tip and the sample. 
When graphene sheets [5] or TI surfaces [5] are stud- 
ied by STM spectroscopy, the experimental DOS curves 
are never simple F-shaped curves. This is unsurprising 
since disorder, electron-phonon, and electron-electron in- 
teractions [7j can all play a role in altering the tunneling 
DOS. Brar et al. have however recently discovered 
some very specific gate- voltage dependent, i.e. carrier 
density dependent, features in the rf/tip/dVbias spectra 
of doped graphene sheets deposited on Si02 which they 
attributed to electron-plasmon interactions. Similar fea- 
tures with an enhanced amplitude have been observed 
by the same group for graphene sheets deposited on h- 
BN [5] , suggesting that the effect is rather universal. 

In this Rapid Communication we propose electron- 
electron interactions as the source of these DOS spectral 
features. Indeed theory has already predicted [101 HI] 
that electron-plasmon interactions have a large impact 
on the one-body spectral function A{k, io) of doped 



graphene. After carrying out extensive studies of 
A{k,uj) in quasi-freestanding doped graphene sheets on 
(hydrogen-terminated) SiC using angle-resolved photoe- 
mission spectroscopy (ARPES), Bostwick et al. [12] re- 
cently discovered experimentally that coupling between 
electrons and plasmons leads to a substantial reconstruc- 
tion of the MDF conical spectrum. The reconstruction is 
related to the appearance of new composite quasiparti- 
cles, known as plasmarons |13j . which are expected to ap- 
pear when interactions between electrons and the collec- 
tive charge-density (plasmon) oscillations of an electron 
liquid |14j are strong enough. Plasmarons are signaled 
theoretically by multiple solutions of the one-particle 
Dyson equation, and experimentally by satellite bands 
in the ARPES spectra. Since the tunneling DOS is given 
by an integral over all momenta of A{k,uj), plasmaron- 
related features must also be present in the STM spectra. 

Fig. [l^) summarizes our theoretical results for the tun- 
neling DOS of electrons with MDF bands which were ob- 
tained using a continuum model GqW (random phase) 
approximation [TOl E] , with the dimensionless Coulomb 
coupling constant aco — /(hve) = 0.5 and the flavor 
number A^f = 4 set to values appropriate for graphene 
on h-BN. The flavor number accounts for graphene's 
spin and valley degeneracies. {Nf ~ 1 for TI sur- 
faces.) The DOS differs from the non-interacting result 
v^{E) — Ni\E\/{2iTh?v'^), in three ways: i) at large en- 
ergies the DOS is linear, but its slope is altered because 
the quasiparticle velocity v* is renormalized |15l 116] ; ii) 
the DOS does not vanish at any energy; and iii) the 
DOS displays a rich structure of plasmaron peaks sep- 
arated by local minima located near the energies of the 
Dirac crossings 12 in the corresponding ARPES spec- 
trum [Fig.[ljD)]. 

In what follows we introduce and discuss qualitative 
features of the continuum model Hamiltonian and the 
random phase {GqW) approximation that we have em- 
ployed in combination to approximate 2D MDF interac- 
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FIG. 1: (Color online) Panel a) Tunneling density-of-states 
(DOS) Vcciui) of electrons with a MDF band Hamiltonian as 
a function of energy u) measured from the chemical potential 
fi. The DOS is nearly independent of carrier density when it 
is expressed in units of the non-interacting DOS at the Fermi 
energy fc'o(eF), and lj is expressed in units of ep. These results 
are for Ni = 4, n = 10^^ cm~^, and aee = 0.5. Panel b) A 
2D color plot of the spectral function A{k, cj) for the same 
parameters used in panel a). The energies -Eo (long-dashed 
line), El (short-dashed line), and E2 (dash-dotted line) define 
a diamond-like shape in (fc, uu) space [T^ in which the influence 
of electron-plasmon coupling is strongest. The DOS has mea- 
surable spectral structure in the energy range E2 < lj < Eq. 



tion physics to obtain the results illustrated in Fig. [T] We 
then compare the roles of electron-electron and electron- 
disorder interaction in filling in the Dirac model DOS V, 
concluding that interaction effects play the dominant role 
in high quality graphene samples. 

The Coulomb MDF model: We describe 2D electron sys- 
tems with MDF bands using the following Hamiltonian 



n 



(1) 



Here v is the bare Fermi velocity , S is the sample area, 
and ■0^ {''Pk,a) creates (destroys) an electron with mo- 
mentum k and two- valued (pseudo)spin index a, and 
flavor indices are implicit. In the case of graphene, 
a — A, B ifi 'A sublattice-pseudospin index, whereas in 
the case_of TI surface states it is a true spin index. 



In Eq. 1 = Y.t. 



Jpk^a is the density oper- 



ator, and Vq ~ 27re /(eg) is the 2D Fourier transform of 
Coulomb potential. Here e depends on the dielectric en- 
vironment surrounding the 2D MDF fluid. In the case of 
graphene e — [ei + e2)/2, where ei (£2) is the dielectric 
constant of the medium above (below) the sheet. For the 
surface states of a thick TI e = (ei -I- eTi)/2 where exi 
is the TI bulk dielectric constant. Because chalcogenide 
TPs have large dielectric constants, interaction effects on 
their properties will tend to be weak. For this reason 
we concentrate on graphene MDF's in what follows. In- 
teraction effects can be made stronger in chalcogenide 
TI MDF's by preparing thin film samples, but interac- 
tions between top and bottom surfaces will then play an 
essential role. 

Below we describe an electron-doped system with an 
excess electron density n; the properties of hole-doped 
systems can be obtained by appealing to the model's 
particle hole symmetry as we explain below. The cor- 
responding Fermi wave number and Fermi energy are 
fcp = a/ Airn/Nf and ep = vk-p, respectively. In calcu- 
lating the DOS it will be convenient to use the Fermi 
level rather than the Dirac point as the zero of en- 
ergy so that the Dirac-band single-particle energies are 
^k,x = £fc,A — £f = Aw/c — Ep where A = ±1 distinguishes 
the conduction and valence bands. Note that the band 
and interaction terms scale in the same way with changes 
in length scale. The spectra we calculate will therefore 
be independent of fcp if energies are measured in terms 
of ep, apart from weak logarithmic dependences on a 
TT-band width ultraviolet cutoff scale which must be in- 
cluded in the theory. It follows in particular that the 
DOS at the Dirac point is accurately proportional to the 
Fermi energy. 

GqW approximation spectral functions: The DOS is re- 
lated to the spectral function A{k,uj) by 



.{to) = Nf 



+ CX3 



dE iyoiE)A{E/v,Lu) 



(2) 



where A{k,u!) — J2\=±-^>^i^^to) and (suppressing the 
fc, uj variables for simplicity) 
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The quantity T,\{k,ijj) denotes the electron-electron in- 
teraction contribution to the quasiparticle self-energy. 
The real part of 'E\{k,uj) is measured from its value at 
the Fermi surface, which physically represents the inter- 
action contribution to the chemical potential, /iint- Cor- 
respondingly, the variable uj in Eqs. [and also 
in Eq. Q below] represents energy measured from the 
chemical potential ^ = Ef + /iint of the interacting sys- 
tem and in units of the Fermi energy e-p. Physically, 
A\{k,uj) represents |14j the probability density for in- 
creasing or decreasing the energy of the iV-particle sys- 
tem by w upon adding or removing a single particle in 
state |fc,A). An STM experiment probes both occupied 
(w < 0) and empty (w > 0) states, whereas ARPES 
probes the occupied portion of the spectrum only. 

For the quasiparticle self-energy we use the GqW 
(or random phase) approximation [TIIHT^ fTl] . in 
which the self-energy is expanded up to first order 
in the dynamically-screened interaction Wcciq,^^) — 

Vq/£BPA{q,^)-- 

QmT.xik,uj) = Yl j ^™ [^cc('7,f^fc+,A')] 

(4) 

whereJlfc^A = a;-^fc,A, -^fcA,fc'A' = [1-l-AA' cos((/3fc-<^fc')]/2 
is the MDF model chirality factor, and fe+ is a shorthand 
for k+q. The RPA dielectric function is e-RPAil, w) = 1 — 

VqX'ppiq,!^), where XppiQ,^) is the well-known Lindhard 
response function of non-interacting 2D MDFs \17\ at 
arbitrary doping n. The real part of the self-energy is 
obtained from Eq. Q by a Kramers-Kronig transform. 

Note that, since the MDF model is particle-hole 
symmetric, the real and imaginary [5fte T,\{k,uj) and 
3m Ti\{k,uj)] parts of the quasiparticle self-energy for 
electron doping are equal to — Sfte T,-x(k,—uj) and 
3™ T._xik,~(^) for hole doping. Thus A\{k,Lj) and 
Vcci'-^) for electron doping are equal to A-x{k, —u)) and 
z/ee(— w) for hole doping. 

Interaction Strength and Disorder Dependence: As seen 
in Fig. [T]d), the plasmaron bands loose spectral weight 
as they approach the Fermi energy but are otherwise re- 
markably well defined, providing a lower energy echo of 
the quasiparticle Dirac point. The diamond- like shape in 
{k, uj) space JJJ that can be clearly recognized in Fig.[T]D) 
is bounded by plasmaron and quasiparticle peaks which 
cross at a finite value of k at an energy Ei which lies be- 
tween the quasiparticle band crossing at energy Eq and 
the plasmaron band crossing at energy E2 ■ The diamond 
feature represents a massive reconstruction of the MDF 
spectrum. As shown in Figs. [I^) and [2] it has a strik- 
ing impact on the DOS at low energies. Note that for a 
n-doped system the strongest plasmaron features appear 
in the filled portion of the spectra. The particle-hole 
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FIG. 2: (Color online) Panel a) DOS of a 2D interacting MDF 
system [in units of the non-interacting DOS at the Fermi en- 
ergy z^o(eF)l as a function of energy lo (measured from the 
chemical potential and in units of the Fermi energy £p). The 
results plotted in this figure are for A'f = 4, carrier density 
n — 10^^ cm~^, and for three values of the coupling constant 
Oec (including the non-interacting case Oee = 0). The vertical 
dotted line represents the chemical potential pt of the inter- 
acting system. Panel b) Expanded view of the energy region 
near the reconstructed Dirac point. The features shown here 
originate from (fe,tj) diamond explained in the text and illus- 
trated in Fig. [T]d) . In this panel we have also included data 
for Qcc = 2.2, a value that is suitable to describe a suspended 
graphene sheet. Note that the diamond features in Ucc{i^) 
are strongest for Oec ~ 0.5, a value thought to be appropri- 
ate to a graphene sheet on h-BN, and not for the strongest 
interactions considered. 



symmetry properties explained above imply that the op- 
posite occurs for a p-doped system. At large energies 
we clearly see that the DOS grows linearly with a slope 
which decreases as Ofee increases, an effect which can be 
understood qualitatively in terms of the interaction en- 
hancement fl6' of the quasiparticle Fermi velocity v* . 
In Fig. [2] we plot the DOS of an interacting 2D MDF 
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Born approximation (SCBA) [TS] for which the disorder 
self-energy reads 



FIG. 3: (Color online) DOS of a 2D non-interacting MDF 
system at the Dirac-point energy in the presence of Coulomb 
scatterers with Z = \ located on the 2D plane where elec- 
trons move [d = in Eq. ([5|]. The dimensionless parameter 
on the horizontal axis is the impurity concentration, i.e. nimp 
in units of the electron density n. The DOS disorder back- 
ground VAm = i-'ci(O) is plotted in units of the non-interacting 
MDF clean-system DOS at the Fermi energy, i.e. i'o(eF), to 
facilitate comparison with Fig. [2] The three sets of data refer 
to three values of the coupling constant Qee- The solid lines 
are guides to the eye. The thick (blue) crosses indicate the 
values of impurity concentration at which j/dis equals i/bg, the 
background DOS of the clean interacting system at the same 
value of Qcc. The two contributions are clearly comparable 
only for impurity concentrations well above unity. The inset 
shows vc\{E) as a function of energy E measured from the 
Dirac point (and in units of the Fermi energy). Results in 
this figure are based on the SCBA - Eq. ([5|. The abruptness 
of the DOS increase is certainly an artifact of the SCBA. 



system for three values of the coupling constant aoc- 
Due to the increase of screening with increasing dee; the 
amplitude of the plasmaron peaks in the energy range 
E2 < oj < Eq becomes quite small when ctco ^ 2. The 
DOS does not vanish at any energy when interactions 
are included, as observed experimentally in Ref. [H The 
strength of the interaction induced background in the 
non-interacting V can be characterized by the absolute 
minimum value of i'oc(w): v^g = mintj[i'cc(w)]. To esti- 
mate how important this background typically is for the 
interpretation of experiment, we can compare it with the 
finite DOS background due to disorder. To this end, we 
calculate the DOS of a 2D MDF system in the presence 
of a model disorder potential, including electron-electron 
interactions only through the static RPA screening func- 
tion eRPA(<;, 0). For the sake of definiteness, we consider 
Coulomb impurities (CI) located at random in a plane 
parallel to the 2D electron system that is separated by 
distance d. Let nimp be the impurity density. For our 
qualitative goals it is sufficient to use the self-consistent 
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, (5) 



where Wci{q,0) = Vei{q)/eYiPA{q,0) is the statically- 
screened electron- impurity potential, and Vei(q,0) = 
VqZ exp{—qd) is the 2D Fourier transform of the bare 
electron-impurity potential (assuming that all impurities 
have charge equal to Ze, in absolute value). The DOS 
in the presence of disorder, vciiE), can be computed by 
iteratively evaluating the integral in Eq. ([5]) numerically 
until T,\ {k,E) is self-consistent. 

In Fig. [3] we plot the DOS background due to long- 
range scatterers, i^dis — i^ci(0)j as a function of the im- 
purity concentration rii^p/n and for different values of 
the fine-structure constant aco- Note, indeed, that since 
we have assumed Coulomb disorder, the strength of the 
electron-impurity interaction is controlled by the same 
dimensionless coupling constant which determines the 
strength of electron-electron interactions. We clearly see 
from this plot that the impurity concentration needed 
to give a background comparable to v^g is rather high. 
For example, for graphene on Si02 the impurity den- 
sity riimp necessary to make disorder more important 
than interactions is more than twice the electron density. 
Typical impurity densities, extracted from a comparison 
between theory and transport data |19) which assumes 
Coulomb scatterers, are for comparison on the order of 

10^° — 10^^ cm~^, giving typical values of rUmp/n < 1. 
We thus conclude that the background due to electron- 
electron interactions is very significant, at least in com- 
parison with that due to long-range Coulomb impurities. 
Of course, the disorder effect will always dominate at 
very small carrier densities. 

In summary, we have calculated the tunneling DOS of 
a 2D system of interacting MDFs. We have found that 
plasmarons in the one-body spectral function give rise to 
a non-trivial structure of peaks and valleys in the DOS. 
The spectral weight carried by plasmarons gives rise to 
a finite background in the DOS, which does not vanish 
at any energy and instead displays two peaks associated 
with the electron and plasmaron Dirac points [12) . A 
more direct comparison with the experimental results in 
Ref. O is necessary to determine convincingly whether 
or not the DOS peak has already been observed in the 
STM spectra of graphene on h-BN. Our predictions can 
also be tested by low-temperature STM spectroscopy on 
quasi-freestanding doped graphene sheets on (hydrogen- 
terminated) SiC [5D]. Finally, our findings are also rel- 
evant to STM studies of TI surface states. A finite 
background has, for example, recently been observed by 
Alpichshev et al. [5] in the STM spectra of the surface 
states of Bi2Te3 which could be due in part to interac- 
tions. 
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